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COMPARISON OF SINK PARTICLES IN AMR AND SPH 

Christoph Federrath^'^, Robi Banerjee\ Paul C. Clark\ and Half S. Klessen^ 

draft February 28, 2010 

ABSTRACT 

Star formation is such a complex process that accurate numerical tools are needed to quantitatively 
examine the mass distribution and accretion of fragments in collapsing, turbulent, magnetized gas 
clouds. To enable a numerical treatment of this regime, we implemented sink particles in the adaptive 
mesh refinement (amr) hydrodynamics code flash. Sink particles are created in regions of local 
gravitational collapse, and their trajectories and accretion can be followed over many dynamical 
times. We perform a series of tests including the time integration of circular and elliptical orbits, the 
collapse of a Bonnor-Ebert sphere and a rotating, fragmenting cloud core. We compare the collapse 
of a highly unstable singular isothermal sphere to the theory by Shu (1977), and show that the sink 
particle accretion rate is in excellent agreement with the theoretical prediction. 

To model eccentric orbits and close encounters of sink particles accurately, we show that a very 
small timestep is often required, for which we implemented subcycling of the A^-body system. We 
emphasize that a sole density threshold for sink particle creation is insufficient in supersonic flows, if 
the density threshold is below the opacity limit. In that case, the density can exceed the threshold 
in strong shocks that do not necessarily lead to local collapse. Additional checks for bound state, 
gravitational potential minimum, Jeans instability and converging flows are absolutely necessary for 
a meaningful creation of sink particles. 

We apply our new sink particle module for flash to the formation of a stellar cluster, and compare to 
a smoothed particle hydrodynamics (sph) code with sink particles. Our comparison shows encouraging 
agreement of gas properties, indicated by column density distributions and radial profiles, and of sink 
particle formation times and positions. We find excellent agreement in the number of sink particles 
formed, and in their accretion and mass distributions. 

Subject headings: accretion, hydrodynamics, ISM: kinematics and dynamics, methods: numerical, 
shock waves, stars: formation 



1. INTRODUCTION 

Molecular clouds are turbulent, magnetized, self- 
gravitating objects. Supersonic turbulence, in par- 
ticular, plays an important role in shaping the cloud 
structure, and in controlling star formation, because 
it creates the seeds for local gravitational c o llapse 



jElmcEreen & Scalo' '20p; 'Scalo & Elmegreen' '20M; 
[Mac Low & Klcsscn 2004; McKcc & Ostriker 2007). 
Due to the filamentary, fractal structure of the interstel- 
lar medium, and due to the large density contrasts, star 
formation in turbulent molecular clouds proceeds in mul- 
tiple regions at the same time in parallel, reflecting the 
hierarchical and fractal nature of th e gas density proba- 
bility d istribution (e.g., iScalol 119 90': 'Vazquez-Semadeni' 
1994 ; lElmegreen fc Falgarond 1996: Klessen et al.i 



To model this complex interplay of turbulence and 
gravity that eventually leads to cloud fragmentation and 
to stellar birth with a well-defined initial mass function, 
it is necessary to follow the freefall collapse of each in- 
dividual fragment, while keeping track of the global evo- 
lution of the entire cloud at the same time. The funda- 
mental numerical difficulty with this approach is that the 
freefall timescale iff decreases with increasing density, 
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stars typically 
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form in clusters (jLada fc Ladal I2003D . showing similar 
fractal patterns as the gas clouds from which they form 
(jSanchez fc Alfaro 2009) . 
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Following the freefall collapse from typical molecular 
cloud densities up to stellar densities requires the nu- 
merical scheme to cover about ten orders of magnitude in 
timescales. Modeling each individual collapse over such 
an enormous dynamic range, and following the large- 
scale evolution over several global freefall times in a single 
hydrodynamical simulation is beyond the capabilities of 
modern numerical schemes and supercomputers. Thus, 
if one wants to model the large-scale evolution of the 
molecular cloud alongside the collapse of individual re- 
gions far beyond the collapse of the first object, the in- 
dividual runaway collapse must be cut-off in a controlled 
way, and replaced by a subgrid model. 

There are two different subgrid models to tackle this 
problem in numerical simulations. The first approach is 
to heat up the gas that exceeds a given density thresh- 
old, pres, which is necessarily related to the achievable 
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numerical resolution. We call this procedure 'Jeans heat- 
ing'. This relative heating of the gas is often modeled by 
changing the effective equation of state above the den- 
sity threshold (e.g., by setting the adiabatic exponent to 
7 > 4/3 for p > Pics)- Heating up the gas above a den- 
sity threshold increases the sound speed locally and thus 
increases the Jeans length. 
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until the gas is stabilized against gravitational collapse. 
The problem with the Jeans heating approach is that any 
parcel of gas above a given density threshold is heated ar- 
tificially, although the actual gas equation of state should 
still be close to isothermal so long as the density thresh- 
old is below the opacity limit. This is the case for den- 
sity threshol ds smaller than about 10~ ^ -^gcm~'^ (e.g., 
Larson 196^; lPenstonlll969l : iLarsonl 120051 : iJappsen et al.l 
2005i . and references therein). Moreover, gas can become 



denser than the threshold value in shocks that do not nec- 
essarily lead to the formation of a gravitationally bound 
structure. Thus, shocked gas not going into freefall col- 
lapse will be heated up artificially. An additional prob- 
lem of the Jeans heating approach is that the increasing 
sound speed in the h eated regions can reduce the Courant 
timestep (see § 12. 5p to prohibitively small values. 

The second type of subgrid model is to use 
so called 'sink parti cl es', a method invented by 
iBate. Bonnell. fc Pried (|1995D for Smoothed Par- 
ticle Hydrodynamics (sph), and first adopted for 
Eulerian, Adaptive Mesh Refi nement (amr) by 
iKrumholz. McKee. fc Kleh] (|2004D . If the gas has 
reached a given density, a sink particle is introduced, 
which can accrete the gas exceeding the threshold, 
without altering the thermal physics. However, sink 
particles are supposed to represent bound objects that 
are/or will be going into freefall collapse, and thus, a 
density threshold for their creation is insufficient. As for 
the 'Jeans heating', shock compression can temporarily 
create local densities larger than the threshold without 
triggering gravitational collapse in that region. Previous 
grid-based implementations of sink particles are mostly 
based on a density threshold criterion. If the density 
threshold for sink particle creation is smaller than 
the opacity limit (about 10~^^gcm~^), we show that 
spurious sink particles are created in shocks that did not 
create a gravitationally bound and collapsing structure. 
Here, we present an implementation of sink particles for 
the Eulerian, AMR code FLASH that avoids this problem 
by using a series o f checks for gravitational collapse sim- 
ilar to Ba te et al.l (fl995) prior to sink particle creation, 
such that only gravitationally bound and collapsing gas 
is turned into sink particles. 

We show that the star formation efficiency and 
the number of fragments is overestimated, if addi- 
tional, physical checks (e.g., checks for a gravitation- 
ally bound and collapsing state) in addition to a den- 
sity threshold are ignored prior to sink particle for- 
mation. We believe that it is also crucial to inves- 
tigate the resulting mass distribution of the sink par- 
ticles, because any successful numerical and analyti- 
cal model of cloud collapse and star formation is ex- 
pected to account for the observed mass distribution of 



cores and stars, i.e., the clump, co re and stellar initial 
mass functions (e.g.. iKlessen et al][l 998; Klessen 200i; 
lKroupa"2001'; 'Pad oan fc Nor dlund'^2 002c iChabrien | 2003 : 
Mac Low & KlessenI 120041: Elmcgr eeii fc Scald I2OO4I: 



Larson 2005: Bonnell fc Bate 2006: McKee & Ostriker 
1^07; Krumholz fc Bonnell 2007; Hennebelle & Chabrie" 
I2OO8, 2009). The sink particle implementation presented 
here enables us to address the star formation efficiency 
and rate, as well as the mass distribution of fragments 
obtained in numerical experiments in a robust and quan- 
titative way. 

We test our sink particle implementation against a 
number of standard fragmentation and orbit integration 
tests. Since stars typically form in dense cluster, close 
encounters of stars are common, and thus we dedicate 
significant attention to testing our scheme for its abil- 
ity to capture close orbits and close encounters. More- 
over, in the densest existi ng star clusters, even mer ging of 
stars might be possible (jZinnecker fc Yorkil2007[ ). Our 
scheme also supports sink particle merging for such ex- 
treme cases. 

We compare our sink particles with the original sink 
particle implementation of Bate et al. (119951 ) in their sph 
code. It is a standard SPH code, and most SPH implemen- 
tations of sink partic les are based on their approach (e.g., 
IJappsen et al.ll2005l) . We find that flash and sph show 
encouraging agreement in the obtained mass distribution 
of sink particles. This code comparison strengthens our 
confidence in numerical calculations of collapse and frag- 
mentation, using AMR on the one hand and SPH on the 
other hand. Here, we primarily introduce the new sink 
particle module for flash and present a series of initial 
tests for follow-up studies using sink particles. For in- 
stance, there are various SPH studies with sink particles of 
purely hydrodynamic collapse. Our comparison with sph 
shows that these results are robust. However, one would 
also like to test the infiuence of magnetic fields and am- 
bipolar diffusion on star formation as well. The first sig- 
nificant steps toward this were taken in sph simulations 
of m agnetized turbulent clouds recently (Price fc Bati 
l2008f l. However, the numerical representation of tangled 
and wound-up magneti c fields in sph turbulence simu - 
lations is (still) limited (|Brandenburgjl2010HPricell2010D . 
Our grid-based approach of sink particles should allow 
a refined modeling of collapsing, turbulent, magnetized 
clouds in follow-up studies, including additional physical 
processes like ambipolar diffusion and self-consistent jet 
and outflow formation. 

In §[2]the new sink particle implementation for the AMR 
code FLASH is explained in detail. A series of simple and 
more complex tests of the sink particle implementation 
is presented in § [3l In § [4] we analyze the formation 
of a star cluster, and compare column density images, 
radial profiles, and mass distributions of sink particles 
obtained with the FLASH code and with a standard SPH 
code. In § O we summarize our results and discuss the 
importance of checks for gravitational instability in ad- 
dition to a density threshold to avoid spurious creation 
of sink particles, and to avoid overestimating the star 
formation efficiency and rate. 

2. NUMERICAL IMPLEMENTATION OF SINK 
PARTICLES IN THE FLASH CODE 
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2.1. The Basic flash Code 

The FLASH codeEl ()Frvxell et al.l [2OOOI: IDubev et al.1 

is an adaptive mesh refinement (amr) code 
(|Berger fc Colellal 119890 . For purely hydrodynamic 
studies, it uses the pi ecewise parabolic method 
(jColella fc Woodwardlll984f) by default to integrate the 
equations of hydrodynamics. For magnetohydrodynam- 
ical (MHD) studies, flash provides an 8-wave Roe 
solver. In addit ion, a new approxini a te Ri emann solver 
for ideal-MHD (jBouchut et al.l l2007l [2009I ). which pre- 
serves positive states in highly supersonic M HD tur- 
bulenc e was recently developed for flash by iWaaganI 
(j2009( ). The corresponding scheme for preserving pos- 
itive states in purely hydrodynamic st u dies h as been 
teste d successfully in iKlingenberg et al.l (j2007l ). More- 
over, iDufhn fc Pudrita l 2008| ) have recentlv developed 
a non-ideal MHD scheme to model ambipolar diffusion. 
This module is also implemented in the flash code 
and works (like all other physics modules) within the 
AMR framework, flash is parallelized with mpi, and 
output files are written in the versatile hdf5 format. 
The self-gravity of the gas is treated with an iterative 
multigrid solver (here we use the multigrid solver imple- 
ment ed in flash v2.5, recently refined for v3 bv iRickerl 
|2008() . Moreover, a tree-based gravity solver was devel- 
oped for flash (Richard Wiinsch 2009, priv. comm.), 
which is currently being modified to run on graphics 
processing units. Our sink particle implementation is 
compatible with existing flash modules, i.e., it can be 
used with the different hydrodynamical and magnetohy- 
drodynamical schemes, including the ambipolar diffusion 
module, and with either the multigrid or the tree gravity 
solver. FLASH has b een extensively tes ted against labo- 
ratory experiments (Calder et al.ll2002f) and other codes 



integer indexes (i,j, fc), such that 

ijk 
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(iDinionte et al.ll20o4:lHeitmann et al.ll2005l:IAgertz et al.l 



l2007HTasker et al.l 120081: iKitsionas et al.ll2009[ ). 

For the implementation of the new sink particle mod- 
ule, we made use of the iV-body capabilities of flash. 
Once created, sink particles are free to move within the 
Cartesian computational domain, independent of the un- 
derlying grid, i.e., they move in the Lagrangian frame of 
reference, while the grid points are fixed in space (Eule- 
rian frame of reference). The outer boundary conditions, 
i.e., outflow, reflecting, or periodic also apply to the sink 
particles in the simulation box. 



2.2. Sink Particle Creation 

Prior to sink particle creation, it is necessary to per- 
form a number of tests, since we want to avoid creating 
spurious sink particles in regions that are not undergo- 
ing freefall collapse. The basic idea is to first check each 
computational cell for whether it exceeds a given density 
threshold pios- If this is the case, a roughly spherical 
region with a given radius race centered on that cell is 
temporarily created from the gas. This radius is usually 
the same as the accretion radius of the sink particle, so 
we will also call it race in the following (see § 12.31 for the 
implementation of accretion) . We denote the region sur- 
rounding the cell with p > p^^s the control volume V. It 
is defined such that it covers all computational cells with 

^ http : //flash . uchicago . edu/website/home/ 



for aU (iAx)"^ + [jl^yf + (kAz)"^ < rl^^. The central ceU 
of each temporary control volume is at {i,j,k) — (0,0,0), 
and AV{i,j,k) — AxAyAz is the computational cell 
volume at spatial position {i,j,k). It is then checked 
whether the gas in the control volume V 

• is on the highest level of refinement, 

• is converging, 

• has a central gravitational potential minimum, 

• is Jeans-unstable, 

• is bound, 

• is not within race of an existing sink particle. 

These checks ar e similar to the checks introduced in 
iBate et aD (|1995l ) prior to sink particle formation in SPH. 
Only after these conditions are fulfilled altogether, a sink 
particle is formed from the gas within the control volume 
and placed in the center of mass of the gas from which 
the particle forms. Each of the criteria for sink parti- 
cle creation is discussed at more detail in the following. 
The order in which these checks are performed does not 
matter, however, during code development and tests it 
became clear that it is useful to do the least computa- 
tionally expensive checks first to make a preselection of 
cells prior to the more expensive checks. 

2.2.1. Density Threshold 

We introduced sink particles in the flash code to 
follow collapse calculations for many dynamical times 
with out violating the Tru elove criterion for the gas den- 
sity ()Truelove et al.llT997l ). The Truelove criterion states 
that in order to avoid spurious fragmentation in numer- 
ical collapse calculations in grid codes, the Jeans length 
(eq. [2]) mus t be resolved with at least four grid ceUs, 
A.t/Acc > 4 (iTruelove et al.lll997| ). For mhd calculations, 
iHeitsch et al.l (|200lD find that more than four cells are 
required . It should be eniphas ized t hat the resolut i on cri - 
teria by ITruelove etall (|1997[) and IHeitsch et all ()200lD 
are only meaningful in regions that are undergoing self- 
gravitational collapse (further discussed in § I2.2.8|) . Usu- 
ally, adaptive mesh refinement (amr) is used to guar- 
antee that the resolution is always sufficient to satisfy 
these criteria. However, for collapse calculations involv- 
ing multiple collapsing regions the pure amr approach 
only works for the first object going into collapse. This 
is because freefall collapse is a runaway process in which 
the first gravitationally bound over-density collapses the 
fastest. Due to the Courant condition (see 12. 5p . this 
leads to smaller and smaller timesteps, which stalls the 
evolution of the entire simulation. Introducing sink par- 
ticles in regions that are going into freefall collapse on the 
highest level of refinement provides a way for cutting-off 
this runaway process in a controlled fashion. However, 
as discussed in the introduction, it is insufficient to form 
sink particles solely based on a density threshold in su- 
personically turbulent gas. Additional checks are neces- 
sary. 
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2.2.2. Refinement Check 

The Jeans refinement criterion discussed in the previ- 
ous subsection is also used to resolve the Jeans length of 
the gas up to the highest level of the AMR grid hierarchy. 
Only when the Jeans refinement reaches the highest AMR 
level, sink particles are allowed to form. Since the accre- 
tion radius race of a sink parti cle is coupled to the grid 
resolution criterion (see § 14.2.11) , sink particles should al- 
ways be located on the highest level of the AMR hierarchy. 
This is taken care of by an additional refinement crite- 
rion for sink particles, which guarantees that any grid 
cell located inside the accretion radius of any existing 
sink particle will always be adaptively refined up to the 
highest AMR level. 

2.2.3. Converging Flow Check 

In order to guarantee that the gas supposed to form 
a sink particle is in freefall collapse, we introduced a 
check for convergence of the flow towa rd the center of 
the control volume, V • v < 0. Unlike iKrumholz et al.l 
([2004' ) we implemented this criterion such that not just 
the total divergence toward the central cell must be neg- 
ative, but also that the flows along each of the principal 
axes must be directed toward the center. The converg- 
ing flow check alone is insufficient, because V • v < 
can also be fulfilled by a localized collision of multiple 
shocks that do not necessarily produce a gravitationally 
bound structure. Thus, here we use the converging flow 
check primarily as a preselection of cells that are con- 
sidered for sink particle formation, before other, more 
computationally expensive checks are performed below. 
Only the combination of the converging flow check with 
the gravitational potential minimum, bound state and 
Jeans-instability checks actually guarantees that the gas 
is in freefall collapse. 

2.2.4. Gravitational Potential Mimmum Check 

We guarantee that sink particles can only be created 
if the central cell of the control volume V defined in 
equation ([3]) is the minimum of the local gravitational 
potential 4> inside the control volume. The central cell 
(i, k) — (0, 0, 0) must fulfill the constraint 

0(0,0,0) = min[(/.(i,j,fc)] (4) 

ijk 

for sink particle creation. 

2.2.5. Jeans Instability Check 

If the converging flow check and the potential minimum 
criterion are fulfilled, the control volume is checked for 
Jeans-instability. The thermal energy Eth and the grav- 
itational energy Eg^av of the gas in the control volume 
are calculated as follows: 

^th = ^E^^(*'^'^)^'(»,J,fc) (5) 

ijk 

i^grav=^M(z,J,fc)0(z,J,fc), (6) 

ijk 

where Cs{i,j, k) is the sound speed, 4'{i,j^ k) is the gravi- 
tational potential due to the gas mass inside the control 
volume and 

M(z,j,fc)=p(z,j,fc)Al/(z,j,ft) (7) 



is the mass inside each cell (i, j, k). The relation 

l^gravl > 2£^th (8) 

must hold for sink particle creation, which means that 
the gas exceeds the Jeans mass within the control vol- 
ume. In the case of magnetic fields this criterion is modi- 
fied such that the magnetic pressure is taken into account 
as well. The magnetic energy, 

£;mag= ^ 5]|B(z,J,fc)|' Ay(i,J,fc), (9) 
OTT ^ — ' 

ijk 

is computed and added to the right hand side of equa- 
tion ([8]) to get a modified version of the Jeans criterion, 
which takes into account the additional pressure provided 
by magnetic field fluctuations. 

2.2.6. Check for Bound State 

For successful sink particle creation, the total gas en- 
ergy inside the control volume must be negative, 

^'grav + -Eth + i?kin + £^mag < . (10) 

The gravitational, thermal, and magnetic energies are 
computed from equations ©, © and ©. The kinetic 
energy 

ijk 

is determined from the velocity dispersion of the gas, 
where the center of mass motion, 

Y.ijkM{i,j,k)^f{i,j, k) 
Y.l,kM{^.3.k) 

is subtracted. Only if equation (jlOp holds, the gas within 
the control volume is a bound system. 

2.2.7. Proximity Check 

A new sink particle cannot be created within the ac- 
cretion radius of an already existing sink particle. Gas 
that exceeds the density threshold within the accretion 
radius of an existing particle is accreted by that particle, 
if the gas passes the accretion checks as explained in the 
following section. 

2.2.8. Comparison with Other Sink Particle 
Implementations 

Except for the refinement and potential minimum 
checks, all checks used here are analogous to the sink par- 
ticle implementation by Bate et al. (1995) for SPH. Test 
simulations showed that the potential minimum check 
prevents almost all spurious sink particles (particles not 
representing bound and collapsing gas) from forming. 
However, the Jeans instability and bound state checks 
are still required in regions of large velocity dispersions, 
occurring in supersonic, shock-dominated flows. More- 
over, magnetic effects are not taken into account by the 
potential minimum check, but are cover ed with the Jeans 
instability and bound state checks (as in lBate et al.lll995l 
except for our inclusion of the magnetic energy). 

The Eulerian, AMR implementations of sink parti- 
cles by IKrumholz etall (,2004) and ,Wang et al.l pOig) . 
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and the Eulerian, unifor m grid implementation by 
iPadoan fc Nordlundl ()2009D use the density threshold 
check, the refinem ent check fonlv iKrumholz et al.ll206^ 
iWang et aUlMol ) and the converging flow check (only 
the total divergence is required to be negative, which 
is implicitly covered by the density threshold check), 
but no check for potential minimum, Jeans instability, 
bound state and proximity of existing sink particles. 
Thus, a significantly larger number of sink particles is 
typically created , which requires subseque nt me rging in 
iKrumholz eTall (j2004) and Wang ct al. ()2010( ). This 
merging of sink particles is used as an effective accre- 
tion i n addition to their Bondy- Hoyle accretion model, 
while IPadoan fc Nordlundl (f2009l) use no merging of sink 
particles, but direct accretion from the grid based on the 
density threshold. 

ft is important to note here that our si nk particle 
imple mentation (as the implementation by iBate et al.l 
|1995[) allows t he density to exceed t he density thresh- 
old set by the iTruelove et all ()1997f ) criterion in some 
cases. However, this does not mean that our implemen- 
tation violates the Truelove criterion. By definition, the 
Truelove criterion only applies to gravitationally bound 
gas in freefall collapse, for which the density reached the 
threshold due to gravitational collapse, but not due to 
purely hydrodynamical compression alone, typically oc- 
curring in supersonic turbulent flows. In such supersonic 
flows, the gas can reach and exceed the density thresh- 
old in strong shocks (e.g., in isothermal gas the post- 
shock density increases proportional to the square of the 
Mach number). Thus, supersonic turbulent density fluc- 
tuations can cause the gas to exceed the density thresh- 
old. Some of those shocks are dominated by the velocity 
dispersion and not by self-gravity. In such cases it is 
important to avoid spurious creation of sink particles, 
because they would not represent gravitationally bound 
and collapsing objects. Our sink particle implementation 
avoids spurious creation of sink particles by explicitly 
testing the gas for local gravitational collapse. It should 
be noted that the additional checks to the density thresh- 
old are only necessary, if sink particles are created at 
densities below the opacity limit (about 10~^^gcm~^), 
when the gas is still in the roughly isothermal regime, 
as for in stance in the la rge-scale simulations of colliding 
flows bv lBaneriee et all (2009) and i n the t urbulent box 
calculations bv IPadoan fc Nordlundl (|2009f ). where sink 
particles are formed at much lower densities to represent 
clusters of stars rather than individual stars. 

2.3. Gas Accretion 

As soon as a sink particle was created, it can gain fur- 
ther mass by accreting gas from the grid. The combined 
gravitational attraction of the sink particle and the gas 
typically lead to an increase of the gas density over the 
density threshold prcs within the accretion radius race of 
an existing particle. If a cell (i, j, fc) within race exceeds 
Pres, the mass increment 

AM = [p{i,j, fc) - prcs] AV{t,j, fc) (13) 

is calculated. If AM is bound to the collective mass of 
the central sink particle and the remaining gas within 
race, AAf is accreted by that particle. To verify that 
AM is bound to the particle, the kinetic energy of AM 



is calculated in the reference frame of the particle and 
compared to its gravitational binding energy. We fur- 
thermore check that AM is moving toward the sink par- 
ticle, i.e., the radial velocity of AM must be negative. 
This additional check will also allow us to model mass- 
loaded protostellar jets from within the control volume 
in follow-up studies. 

If the mass increment AM is located inside a region 
of overlapping multiple sink accretion radii (which can 
happen due to particle motion) , the gravitational binding 
energy of AM for each of these particles is calculated, 
and AAI is accreted by the particle to which it is most 
strongly bound. 

Gas exceeding the density threshold prcs within a given 
inner accretion radius will always be accreted without 
further checks for a bound state and convergence. The 
inner accretion radius, however, is always set such that 
accretion without further checks only applies to the single 
central cell in which the particle is located. This is only 
to avoid numerical problems like division by zero and 
inflnite gravitational energies during the default checks. 

If the gas mass AM has successfully passed all the 
aforementioned tests for accretion, it is accreted by the 
central particle, such that mass, linear momentum, and 
angular momentum are conserved (see Appendix § [B] for 
a discussion of angular momentum conservation). The 
accreting particle is moved to the center of mass of the 
particle-gas configuration before the accretion step. 

2.4. Gravitational Interactions 

We compute four different contributions to the gravi- 
tational interactions between the gas on the grid and the 
sink particles: 

1. gas-gas (g-g) 

2. gas-sinks (g-s) 

3. sinks-gas (s-g) 

4. sinks-sinks (s-s) 

The modeling of these interactions is described in detail 
in the following sections. In the code, these interactions 
are computed in the order given above, so we also present 
them in this order below. Physically, however, the order 
does not matter. 

2.4.1. Gas-Gas 

The self-gravity of the gas is computed using the stan- 
dard Poisson solver in flash. It is an iterative multigrid 
solver that cycles over the AMR hierarchy to solve Pois- 
son's equation 

V2$gas=47rGpgas (14) 

for the gas distribution on the grid. A tree-based grav- 
ity solver was developed for FLASH recently (Richard 
Wiinsch 2009, priv. comm.), which can also be used to 
calculate the gravitational potential instead of the multi- 
grid solver. 

The Poisson solver returns the gravitational potential 
<l>gas in each grid cell, which is due to the whole gas dis- 
tribution only (without the sink particle contribution). 
The gravitational acceleration of the gas, gg-g is then 
computed for each grid cell as 

gg-g - -V<f>gas . (15) 
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2.4.2. Gas-Sinks 

Equation is used to calculate the gravitational ac- 
celeration for the sink particles due to the gas compo- 
nent only. The acceleration gg-g is interpolated from 
the grid onto the sink particles with a first-order cloud- 
in-cell method at each position of a sink particle n, which 
yields gg-s, n ■ A higher-order interpolation scheme like 
the triangular-shaped cloud method or the tri-cubic in- 
terpolation (e.g., L ekien fc M arsden 2005) does not yield 
significantly different results, because the gravitational 
acceleration of the gas is relatively smooth, and a linear 
interpolation scheme is sufficient. 

2.4.3. Sinks-Gas 

Due to their mass the sink particles can exert an ap- 
preciable gravitational acceleration onto the gas as well. 
This acceleration is computed by a direct sum involving 
all computational cells and all sink particles. For each 
computational cell center (i, j, k), the distance r„ to each 
particle n is calculated and the acceleration that the par- 
ticles with masses M„ exert onto the gas in each cell is 
then computed as 



gs-g(j,J, k) 



E 



GM„ 



(16) 



Computing the acceleration of the gas due to the sinks 
thus involves a nested sum over all grid cells and all par- 
ticles. This can become computationally expensive, be- 
cause this operation scales as the number of grid cells 
times the number of particles. However, even for sink 
particle numbers up to 10'^, the nested sum hardly af- 
fected the overall computational cost, which always re- 
mained dominated by the computation of the self-gravity 
of the gas component (eq. [HI) . 

For cell centers k) very close to sink particles, the 
acceleration computed via equation (I16p can become very 
large and goes to infinity when a sink particle is exactly 
located in the center of a cell. It is thus necessary to 
use gravitational softening of the acceleration within a 
softening radius of each sink particle, which is discussed 
in §[2X3 

The acceleration of the gas due to the sink particles, 
gs-g is added to the gas-gas acceleration gg-g to yield 
the total gravitational acceleration for the gas: 



ihj^k) 



r{hj,k) +gs-g{i,j,k) (17) 



2.4.4. Sinks-Sinks 

The gravitational acceleration for sink particle n due 
to all other sink particles in the domain is computed as 
a direct sum involving all other sink particles m with 
masses M^: 



E 



GM„ 



(18) 



where r^m = r,„ — r„, is the relative distance vector 
between two sink particles n and m. 

The sink-sink acceleration is added to the grid- 
interpolated acceleration caused by the gas (cf. 12.4.21) to 
yield the total gravitational acceleration for sink particle 
n: 

gsinks, n — gg— s, n ~i~ gs— s, n (-^^) 



Equation (|T5)) is also subject to gravitational softening 
as explained in the next section. 

2.4.5. Gravitational Softening 

The basic problem in computing the gravitational ac- 
celeration following equations (fTB|) and (|T8l) is that it 
can yield extremely large values, if the distance between 
sink particles and cell centers or between adjacent sink 
particles becomes small. If the distance goes to zero, 
the acceleration becomes infinite. Hence, the timestep 
given by equation (|23| goes to zero and the simulation 
grinds to a halt. It is therefore necessary to apply gravi- 
tational softening for distances smaller than a given soft- 
ening radius rgoft, such that the acceleration smoothly 
approaches zero as the distance between particles goes 
to zero. There are different approaches for gravitational 
softening. One of the standard softening types is the 



Plummer softening, g(r)cx(|r| 



' soft 



) r. In the case 



of Plummer softening the acceleration approaches New- 
ton's acceleration for a point mass g (x in the limit 
r ^ fsoh ■ Instead of Plummer softening, we use a type 
of spline softening that is typically used in SPH and A^- 
body simula tions to soften the gravitational forces (e.g., 
iPrice fc Mo naghan 2007). The functional form of the cu- 
bic spline softening used here is provided in Appendix [Al 
Figure [T] shows a compar ison of Plummer softening and 
spline softening, equation (jAip . For Plummer softening, 
the gravitational acceleration is modified for all distances 
r. In contrast, the spline softening exactly follows the 
asymptotic solution g oa 1 /r^ for r > Tgoft and smoothly 
approaches zero for r < rgoft- 

2.5. Particle Timestep and Subcycling 

Sink particles are evolved using a variable-timestep 
leapfrog integration scheme. We consider four timestep 
constraints to guarantee a stable numerical solution: 



Atr 



= CcFL min 



Ax 



■i,j,k Vmax(|v(j, j, fc)|, Cs) 



Aigg 


— Qg 




— Cvs 



mm 



At„s = Cgs min 



A 



X 



|ggas(j,i, k)\ 

Ax 



1/2 



2|v„| 

niin(|r„„|, Ax) ^ 

ISsinks, n \ 



(20) 
(21) 
(22) 
(23) 



with CcFL, C'gg, Cvs, Cgs < 1, and Ax is the smallest 
linear cell size in the computational domain. The small- 
est of the first three timestep constraints applies to the 
hydrodynamic solver to guarantee a stable hydrodynam- 
ical evolution of the gas. Equation (l20l) is the Courant- 
Fried richs-Lewy condition ([Courant. Friedrichs. fc Lewvl 
[1921) . A modified version of this is used for magne- 
tohydrodynamic studies, which takes into account the 
fastest possible magnetosonic wave. By solving equa- 
tion ([2T|) , we also consider the gravitational forces of the 
self-gravity of the gas, and the gravitational forces of the 
sink particles exerted onto the gas. Equation (|^^ en- 
sures that sink particles cannot cross more than half a 
grid cell {Ax/2) within a single timestep when they move 
with a velocity v„. The latter two constraints guarantee 
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that the gravitational forces exerted by the sink particles 
onto the gas are consistently taken into account for the 
evolution of the hydrodynamics. 

Since the sink particles are always moved within a hy- 
drodynamical timestep, the sink particle timestep can 
never become larger than the smallest of the hydrody- 
namical timesteps given by equations (1^ . (pij) and (1^ . 
However, these three constraints are insufficient to guar- 
antee a stable and accurate time integration of sink par- 
ticles, because the sink particle timestep can become 
significantly smaller than the hydrodynamical timestep 
for close encounters of sink particles. One additional 
timestep constraint that captures closely approaching 
sink particles is needed. 

The last timestep constraint given by equation (j23p 
takes into account the total gravitational acceleration, 
lgsinks.nl outo cach siuk particle (eq. I19|). and the mini- 
mum distance |r„m| between all sink particles n and m. 
In the course of code development and initial tests of the 
scheme, it turned out that fulfilling equation (|23l) is ex- 
tremely important to avoid inaccurate time integration 
and artificial acceleration, especially in cases of closely 
approaching or orbiting sink particles. In test simula- 
tions and applications, Atgg became up to three orders 
of magnitude smaller than the timestep given by any of 
the three constraints of equations (PH)) . (PTI) and (f^ . To 
avoid prohibitive small timesteps for the hydrodynamical 
evolution, we use subcycling for the time integration of 
sink particles, while keepi ng the hydro dvnamical gas dis- 
tribution fixed fsimilar to lKrumholz et al.. 200 4). During 
subcycling, sink particles cannot change their positions 
by more than half a grid cell due to equation (f^ . and 
thus they remain almost fixed relative to the grid. How- 
ever, subcycling is absolutely necessary to guarantee ac- 
curate integration of sink particle orbits as shown in a 
series of A^-body tests in the next section. 

3. TESTS 

In the following sections we describe a series of test 
simulations to analyze the performance of our sink par- 
ticle scheme. In particular, we test the accuracy of the 
time integration of pure particle systems and systems in- 
volving gas-sink interactions as well as dynamical sink 
particle creation and their accretion. 

3.1. N-body Tests 

As described in I2.4l we modified the gravitational in- 
teractions of gas and particles compared to the standard 
FLASH implementation. Therefore, we test our new ap- 
proach with the following setups. 

3.1.1. Circular Orbits 

In the simplest tests we initialize two sink particles 
with equal masses of 1 M© at a distance of 1 AU, and let 
them orbit around their common center of mass. There 
are 8x8x8 grid cells in this test. However, the gas den- 
sity was set small enough, so that this setup purely tests 
the A^-body integration scheme. The results are shown 
in Figure [2l The trajectory of the particles are plot- 
ted including all positions up to 1000 orbits around their 
common center of mass at (x,?/, z) = (0,0,0). The or- 
bit is maintained circular for these 1000 orbits, and it is 
not expected to deviate significantly for longer integra- 
tion times. This is because the leapfrog time-integration 



scheme is symplectic for constant timesteps. The scheme 
thus conserves energy and preserves time-reversibility in 
this test, such that the symmetry of the system is main- 
tained to machine precision. 

3.1.2. Elliptical Orbits 

A more demanding test is the integration of an ellipti- 
cal orbit. In this case the leapfrog integrator takes vari- 
able timesteps depending on the distance and gravita- 
tional acceleration according to the timestep criterion, 
equation ([23|) . For this test we use the same initial con- 
ditions as for the circular orbit test, but rotate the initial 
velocity vectors by 45° from the tangent of the circular 
orbit of the previous test. The system thus has the same 
total energy, but the sink particles move on elliptical or- 
bits around their common center of mass. The results 
are shown in Figure |3l The top panel shows the default 
time integration with subcycling, which means that equa- 
tions (Uni), (EH), (E21) and (1231) are all used in combina- 
tion to determine the timestep. The bottom panel shows 
the same, but only taking into account the first three 
of these timestep constraints. It is important to note 
that equation (|22|) was the only additional timestep con- 
straint for simulations with particles in the flash code 
before our modifications. Figure [3] (bottom panel) shows 
that using equations (1^ . (PT|) and ([^ is clearly insuf- 
ficient, if close orbits are to be accurately reproduced. 
Even within a single orbit, a significant amount of arti- 
ficial angular momentum is added to the particles if no 
subcycling is performed, and after ten orbits they have 
accumulated an artificial perihel shift of about 21°. In 
contrast, ten orbits are accurately followed without a no- 
ticeable shift with our present scheme. It should be noted 
however that even when using subcycling, the symplectic 
property of the leapfrog integrator is broken, and after 
roughly 100 elliptical orbits, an error leading to a perihel 
shift of about 1° shows up. 

3.1.3. Gravitational Softening and Subcycling Test 

This setup tests the gravitational softening introduced 
in § 12.4.51 We use the same computational domain as 
in the two previous tests with two 1 Mq sink particles 
separated by 1 AU. One particle is initially located at 
{x,y,z) = (—1.0,-1-0.5, 0.0) AU, and the other one is at 
(—1.0, —0.5, 0.0) AU. Both particles have zero initial ve- 
locities in the y-direction, such that they accelerate to- 
ward each other by their mutual gravitational attrac- 
tion. While approaching each other, the timestep drops 
according to the constraints given by equations (|22p 
and (|23p . Due to the softening, however, the timestep 
does not go to zero as the particle distance goes to zero 
(cf. Fig.[T]). This allows the two particles to pass through 
each other smoothly. By the time they pass through 
the point of zero distance, they have converted their ini- 
tial potential energy completely into kinetic energy. The 
softening affects the specific value of this energy but nev- 
ertheless, the energy must be conserved. Thus, after one 
passage both particles have exactly exchanged their ini- 
tial positions. This process then starts again and the 
system oscillates for an infinite time. In order to make 
these oscillations visible, we apply an initial constant 
velocity in x-direction to both particles. The result is 
shown in Figure |4l The top panel shows both particle 
trajectories after 1000 oscillations. Using subcycling, the 
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system easily conserves energy during these 1000 oscil- 
lations. In contrast, if the timestep constraint given by 
equation (j23|) is ignored, energy conservation is broken 
already during the first oscillation as shown in Figure |4] 
(bottom panel). Both particles gain energy artificially 
until they leave the domain. In the current test, this 
happens after about 68 oscillations. Scattering experi- 
ments have shown that violating the timestep condition 
given by equation ((23)) . i.e., not using subcycling can 
lead to enormous artificial accelerations already during a 
single close encounter of two particles. 

3.2. Gas-Sinks Gravity and Refinement Test 

The following setup tests the accuracy of the time in- 
tegration of sink particles orbiting in the gravitational 
potential of a static gas distribution. We create a singu- 
lar isothermal sphere with the following density profile 
on the AMR grid: 



p{r) = p{ro) 



(24) 



,-3 



where ro = 5 x 10 °cm and p{ra) = 3.82 x 10^ ° gem 
The sphere thus has a mass of roughly 3 AIq for r < tq. 
We place three sink particles at radii of r = 1, 2, 3 x 
10^^ cm with masses of 10~^° Mq. The particles thus rep- 
resent test particles in a gravitational potential caused by 
the gas only, such that sink particle interaction is neg- 
ligible. Since the singular isothermal gas distribution 
would immediately collapse and form a new sink particle 
in the center, we deliberately switch off the hydrody- 
namical evolution, and keep the density distribution and 
associated potential artificially static. However, we keep 
the AMR scheme and the Poisson solver fully operational, 
and adaptively refine on Jeans length and sink particles 
as discussed in § 12.2.21 We use a base grid of 32^^ grid 
cells plus two levels of refinement with isolated boundary 
conditions for the gravity. The Jeans length is resolved 
with 6 grid cells for this test, while de-refinement is trig- 
gered as soon as the Jeans length is resolved with more 
than 12 grid cells. The initial density and sink particle 
distributions are shown in Figure [5] (left). 

The three sink particles were given an initial velocity 
in the y-direction, such that they should remain on circu- 
lar orbits around the center of the gas sphere. It is easy 
to determine this velocity analytically from the density 
distribution, equation (j24l) . The result is a Keplerian ve- 
locity of vk = \/GM{r)/r = 0.895 km s~^ for each sink 
particle, independent of the radius. However, since the 
sphere is not singular in the center due to the finite reso- 
lution of the grid, the actual velocities necessary to keep 
the sink particles on circular orbits are slightly smaller 
than the analytical solution (by roughly 1-2%). 

Figured] (right) shows the AMR hierarchy and the sink 
particle positions after the innermost particle has fin- 
ished roughly one orbit. At this time, the sink particle 
on the intermediate orbit with r = 2 x 10^^ cm has fin- 
ished almost half an orbit, and the outer sink particle 
has completed one third of an orbit. The grid is refined 
according to the position of the sink particles. The Jeans 
refinement criterion keeps the central part at the highest 
level of refinement and the outer parts at a resolution 
such that the Jeans length is refined with at least 6 grid 
cells. De-refinement to the base grid resolution was not 



triggered in the outer parts, because the Jeans length 
was resolved with 6-12 grid cells there. 

We evolved this configuration for 100 orbits of the in- 
nermost sink particle. The trajectories of all three sink 
particles are shown in Figure [5] after 1, 10, and 100 orbits 
of the innermost sink particle (from left to right). After 
100 orbits of the innermost sink particle, the intermedi- 
ate particle has finished 50 orbits, while the outermost 
sink particle has finished about 33.3 orbits. The latter 
orbit as well as to a smaller degree also the intermediate 
orbit are slightly broadened, and they are not as accu- 
rately reproduced as the innermost orbit. This is because 
the Poisson solver introduces slight deviations from the 
spherical gravitational potential mainly in regions close 
to the boundaries of the Cartesian domain at the coarse 
grid resolution used in this test. For such a configuration 
and setup that incorporates the multigrid solver and the 
full AMR scheme, it is hard to maintain perfect spherical 
symmetry for more than 10-100 sink particle orbits. 



3.3. 



Collapse of a Bonnor-Ehert Sphere 

The collap se o f supercritic al Bonnor-Ebert 

spheres (lEbertl [l955l: IB onnon [l956h. is a well-studied 
problem (e. g. lLarsonl1l969t iFoster fc Chevalie^ 119931: 



iBaneriee etall 120041: [Aikawa et all I2005D . and thus 
represents a good test case for our sink particle imple- 
mentation. A Bonnor-Ebert sphere can be described by 
the following dimensionless parameters for the radius, 
time, mass, and mass accretion, respectively: 



t 



1/V4FG^ 
M 



m - 



M 



cl/G' 



(25) 
(26) 
(27) 

(28) 



where Cg is the sound speed and po is the central den- 
sity of the sphere. In the hydrostatic configuration of 
such a sphere, the dimensionless mass can be written 
as m = where cj)' — —d\n{p/po)/d^ is the grav- 

itational acceleration. The sphere becomes supercriti- 
cal if its dimensionless radius exceeds the critical value 
Ccrit — 6.451. To trigger the coUapse, we set the initial 
value to ^max = 7.0, which gives a total mass, m = 17.3 in 
dimensionless units. We follow the collapse with two dif- 
ferent density thresholds for sink particle creation and ac- 
cretion: pi-cs = 10^ Po and pros = 10"* po- These two den- 
sity thresholds correspond to accretion radii of ^ — 0.1 
and — 10~^, respectively. 

Figure [7] shows the mass accretion history of the sink 
particle in the case of the low density threshold, pros — 
10^ Po- Initially, when the sink particle starts to accrete 
at Tsink — 0, the accretion rate increases quickly until it 
reaches a maximum. The accretion rate peaks at Tgink ~ 
0.3 when about 12% of the total mass of the sphere is 
accreted onto the sink particle. After that, the accretion 
rate decreases while the envelope runs out of gas. At 
Tsink = 1-62 about 64% of the total gas mass has been 
accreted onto the sink particle. 
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The time evolution of the accretion rate is similar in 
the case of a higher density threshold (pres = 10'* po) for 
the sink particle accretion. Figure [8] shows the mass evo- 
lution for this case. Here, the effective accretion radius is 
smaller by a factor of ten compared to the Pres = 10^ Po 
case, and thus the peak accretion rate is reached ten 
times faster. 

The time evolution of the mass accretion is 
in qualitative agj reemei it wi th ID simulations by 
[Foster fc C hevalier (e.g., '1993"), and with various ana- 
lytical (e.g.. Hunter 1977; Whitworth & Summers 1985) 
and numerical models of protostellar collapse, which pre- 
dict a peak in the accretion fo llowed by a smoothly de- 
clining accretion rate (see also lSchmeia fc Klessen|[200^ 
and references therein). 

3.4. Collapse of a Singular Isothermal Sphere 

The collapse of a singular isothermal sphere with a 
p(r) (X density profile produces a constant flux of 
mass through spherical shells, i.e. a constant accretion 
rate (see e.g. Shu 1977). We compute the collapse of 
a singular isothermal sphere to test whether our model 
of sink particle accretion reflects this collapse behavior. 
The sphere has a truncation radius of i? = 5 x 10^^ cm, a 
density at this radius of p{R) — 3.82 x 10~**gcm~^, 
and therefore a total mass of 3.02 M©. The sphere 
is at a temperature of 10 K corresponding to a sound 
speed of 0.166 kms^^. With these values the sphere 
has a large instabilit y paramet er oi A — 29.3, where 
A = At:Gp{R)R'^/cI IsM [19771) . 

Figure [HI shows the mass and the mass accretion rate 
onto the sink particle during the collapse of the singular 
isothermal sphere. As the gas of the sphere is initially 
at rest, the mass accretion reaches a close to constant 
value of about 1.5 x lO~''M0yr~^, until the entire gas 
of the sphere is accreted onto the central sink particle at 
t~2 X lO^yr. 

For comparison of the accretion history of the sink 
particle, we show the time evolution of radial profiles 
of the gas sphere in Figure [TUl As expected, the den- 



sity profile of the sphere evolves from r " to r " " {see 
e.g.. lShiil[T977tlOgino et al.lll999D . The collapse proceeds 
highly supersonically {M. > 15) with increasing infall 
velocities toward the collapse center (Fig. [TOl middle 
panel). The mass accretion, i.e. the mass flux through 
spherical shells of radius r is constant in time as well 
as independent of the radial position with a value of 
M sa 1.5 X lO~^M0yr~^. It is the exact same amount 
of gas that flows through these spherical shells into the 
control volume of the sink particle, where it is flnally ac- 
creted. A comparison of the sink particle accretion rate 
in Figure [HI with the mass flux through spherical shells 
in Figure [TO] (bottom panel) shows the high accuracy of 
the sink particle accretion mechanism. 

For compar ison, we ca lculated the mass accretion rate 
predicted by IShul (jl977[) . The consta nt accretion rate 
often referred to in many studies citing IShul ()1977[ ) is 



-1.5 



M 



mo 



G 



(29) 



with mo = 0.975, which is the appropriate value for 
an instability parameter close to ^ = 2 (see IShul 119771 
Tab. 1). The accretion rate depends on the coefficient 



mo, which is controlled by A. Numerical integration of 
the dimensionless, one-dimensional equations of hydro- 
dynamics gives the solutions of m o for given instability 
parameters A as in [Shu (see [MS, Tab. 1). There, Shu 
only tabulated values up to A = 4. Since our singular 
isothermal sphere has a much larger instability param- 
eter (A = 29.3), we repeated the numerical analysis by 
IShul ()1977[ ) to find solutions for instability parameters 
up to A = 50. The result for the coefficient mo of the 
Shu accretion rate in equation (|29|) is plotted as a func- 
tion of the instability parameter A in Figure 1111 The 
parameter mo increases from mo — 0.975 for A = 2 to 
mo = 279 for A = 50, more than two orders of magni- 
tude higher! For the instability parameter A = 29.3 of 
the singular isothermal sphere studied here , we obtain 
mo = 133. The accretion rate predicted bv IShul ()1977l ) 
is thus M « 1.45 x 10"''^ M© yr~*, which is in excellent 
agreement with our numerical estimate of the sink par- 
ticle accretion rate (cf. Fig. [Hj). 

3.5. Rotating Cloud Core Fragmentation Test 

We now analyze the collapse and fragmenta- 
tion of a rotating cl oud core, also known as the 
iBoss fc Bodenheimed ()1979| ) test. It is a stan- 
dard test for fragmentation in hydrodynamical codes 
(e.g , iB odcnhcimcr & Boss 1981; Burkcrt ct al. 199]!^ 
iTruelove et al.Hl997l: ICommercon et al.ll2008D. In partic- 
ular w e us e a setup similar to Burkert fc Bodenheimed 
([1991 and IBate fc Burhe^ (119971 1. The initial param- 
eters for the rotating cloud core are: radius i? = 5 x 
10^^ cm, constant density po = 3.82 x lO^*'* gcm^'^, mass 
M = IMq, angular velocity 17 = 7.2 x lO^^^rads"* 
(ratio of rotational to gravitational energy f3 — 0.16), 
sound speed Cs = 0.166 kms~^ (ratio of thermal to grav- 
itational energy a — 0.26), global freefall time tg = 
1.075 X 10^2 g ^ 3 42 X lO^yr, and a 10% density per- 
turbation with an m = 2 mode: p = po[l + 0.1 cos(2</j)], 
where tp is the azimuthal angle. The following polytropic 
equation of state, 

P = c^pr , (30) 
was used with the polytropic exponent 



1 



for p/(10"i5gcm^3) < 0.25^ 



r=0.1 for 0.25 < p/(10-i^gcm-3) < 5.0, (31) 
4/3 for p/(10-i'^gcm-3) > 5.0. 

The initial rotation {/3 = 0.16) forces the gas sphere to 
collapse to a disk with a central bar due to the m = 2 
initial density perturbation. Sink particles are allowed to 
form at densities exceeding 10^*'' gcm^'^. The accretion 
radius was set to race = 39 AU corresponding to 2.5 grid 
cells at the highest level of refinement. 

Figure [T^ shows a face-on column density projection of 
the disk at different times along the collapse of the rotat- 
ing cloud core. Two sink particles form at the location 
of the two fragments at t = 1.26 iff. At t = 1.29 iff a bar 
forms that connects the two main fragments. A third 
fragment forms in the center of the cloud core. After 
that, the two main fragments move toward the central 
object due to the global collapse of the cloud core. At 
t = 1.33 iff two of the three fragments merge to a sin- 
gle particle close to the center of the rotating cloud core. 
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Note that sink particle merging was used in this test. 
Sink particle merging is optional in our implementation 
in FLASH, and can be activated if desired. If activated, 
sink particles are only allowed to merge, if they are inside 
the accretion radii of one another, and if they are grav- 
itationally bound and converging. The merged particle 
is moved to the center of mass of the merging particles, 
and their linear and angular momenta are assigned to 
the merged particle. 

The conservation laws during accretion and merging 
of sink particles are discussed in Appendix |B] The ini- 
tial total angular momentum of the rotating cloud core 
computed from the analytic solution is Lq = 1.44 x 
10^"* gcm^ s~^. The numerical solution conserves angular 
momentum to within 2% for all times including the last 
snapshot shown in Figure [12] when 23% of the mass has 
been accreted onto sink particles. This is comparable to 
the typical conservation of angular momentum achieved 
in numeric al simulations of rotating self-gravitating cloud 
cores (e.g.. lCommer(;;on et al.|[20Q8l ). 



4. STAR CLUSTER FORMATION: AMR VERSUS 

SPH 

To check our sink particle implementation for the AMR 
code FLASH on a more complex physical problem, we ap- 
ply it here to the formation of a star cluster. We further- 
more want to compare our sink particles against an exist- 
ing sink particle implementation in a typica l SPH code. 
We co i npare to th e sph code developed by iBate et aD 
(fT995h . IBate et al.l ^MB) were the first to use sink par- 
ticles in an SPH code, and thus, most SPH implementa- 
tions for sink particl es are based on their approach (e.g., 
iJappsen et al.l [2005( ) . as is the grid-based implementa- 
tion presented here. A detailed description of the FLASH 
code and the our sink particle implementation was given 
in § 12.11 In the following we briefly describe the SPH code 
used for our sink particle comparison. 



4.1. The SPH Code 



Smoothed particle hydrod ynamics (sph, ILucvI 119771 : 
iGingold &: MonaghanI Il977| ) calculations presented in 
this study were performed using a code based on the 
version developed by Benz ()Benj|1988l 119901: iBenz et al.l 
[T990t ). which has since been mo dified by Bate to include 
individual particle timesteps (Nav arro" fc W hite 199^; 
Hern auist & Katz 1989) and sink particles ([Bate et aLi 
19951 ). To capture shocks, the c ode uses the standard 



artificial vi scosity suggested bv [Gingold &: MonaghanI 
(p83) and iMonaghan fc Gingoldi (|1983[ ). with a = 1 
and (3 = 2. To provide adaptive resolution, the code 
allows the SPH smoothing lengths (h) to vary in time 
and space, with the constraint that each particle main- 
tains a roughly constant number of neighbors (50 ± 10) 
within a distance 2h. Gravitational force s are found us- 
ing a binary tree () Press! 1 19871 : lBenz|[l988l and references 
therein), which is also used to obtain the neighbor lists 
required by the SPH algorithm. The binary tree opening 
angle was set to 0.47 in this study, and so returns a more 
accurate calculation of the gravitational forces than the 
theoretical limit of 0.57 for 3D binary trees described by 
I Press! (|1987[ ). Gravitational forces are softened using the 
spline softening technique described in Appendix [X] In 
the SPH code this softening is used between all particle 



types. Between the gas-gas interactions, the softening 
length Tsoft is given by the mean smoothing length of the 
particle pair, and so the softening is zero for particles 
which are not neighbors, since the distance between them 
is greater than 2h by definition. This form of the soft- 
ening assumes that the stand ard SPH smoothing kernel 
(jMonaghan fc Lattanziolfr985D describes the radial den- 
sity profile of the particles, and has been shown to reduce 
the effects of artificial fragmentation (!Bate fc Burkert! 
IT997l:lwifitTOrtl][T99llH^^ et al.ll200fit) . For the sink 
particles, one is free to chose the value of rgoft and in 
our current study this is set to the sink particle accre- 
tion radius, face, as for flash. The value of rgoft be- 
tween gas-sink interactions is then set to (race + h)/2 for 
each pair-wise force. Finally, the SPH equations and sink 
particle trajectories are integrated using a second-order 
Runge-Kutta-Fehlberg integrator. A full descript ion of 
the tim estepping constraints can be found in .Bate et al.l 

dmi). 

4.2. Initial Conditions for the Code Comparison 

An isolated gas sphere of radius R = 5.0 x 10^^ cm = 
0.16 pc with a uniform density of po = 3.85 x 
lO^i^gcm"^ containing M = 100 M© was initialized 
in both the FLASH and the SPH code. The initial 
column density structure is shown in Figure [13] (top 
panels) for the FLASH run (left panels) and the SPH 
run (right panels). An initial random, divergence-free 
turbulent velocity field was generated on a grid with 
128^ grid cells with a velocity power spectrum P{k) cx 
A;"'*, consistent with the observed veloci ty dispersion- 
size relation in mol ecular clouds (e.g., !Larsonl IT98TI : 
!Hever fc Brunti [2004D . and consistent with the velocity 
spectra us ually obtained for driv e n supersonic turbu- 
lence (e.g., [Fed errath etld? '2009a!: !Schmidt et all [2009!: 
!Federrath et al.l .2009^ . However, in the present study 
we are not driving the turbulence, i.e., we study decay- 
ing turbulence. The velocity dispersion was scaled to 
0.89 kms~^, consistent with the observed velocity dis- 
persio ns in molecu l ar clouds of the size studied here 
(e.g., !Larson! !l98l!: !Falgarone erall !1992! ). With the 
isothermal sound speed of Cg = 0.19 kms~^ (tempera- 
ture T = 10 K, mean molecular weight n = 2.3), this 
velocity dispersion corresponds to an initial RMS Mach 
number of about 4.5. Thus, the ratio of kinetic to grav- 
itational energy of the gas sphere is i?kin/-Bgrav ~ 0.25. 
The equation of state is isothermal, and self-gravity is 
used throughout this code comparison. Sink particles 
are allowed to form in both codes at densities exceed- 
ing Pros = 8.0 X 10~^^gcm~'^ with a sink particle ac- 
cretion radius of r^cc = 7.3 x 10^^ cm = 490 AU. Both 
codes used the same spline softening (cf. Appendix [A| 
for sink particle interactions, with the softening radius 
set equal to the accretion radius (rgoft — ^acc)- The tem- 
poral evolution is followed in units of the initial freefall 
time tff = 3.39 x lO^^g ^ i qj ^ 10^ yr. 

For the flash run, we used a cubic computational do- 
main with a boxsize of 0.4 pc, slightly larger than the 
diameter of the initial gas sphere {d = 0.32 pc). To keep 
the mass outside of the sphere small, the gas density was 
set two orders of magnitude smaller outside than the uni- 
form density inside the sphere. The additional mass of 
2.7 M0 outside the sphere is about 3% of the total mass. 
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and thus the total mass in the computational domain is 
3% larger in the FLASH run than in the SPH run. Pres- 
sure equilibrium is provided by keeping the gas outside 
the sphere at a temperature two orders of magnitude 
higher than inside. 

For the SPH run the initial velocity field was inter- 
polated from the grid to the SPH particles using the 
stand ard SPH smoothing kernel (iMonaghan k, Lattanziol 
119851 ). The initial density outside the isolated gas sphere 
is zero for the SPH run since no particles were used there. 
However, a constant boundary pressure was applied in 
order to keep the sphere in pressure equilibrium with the 
surrounding. Both codes used isolated boundary condi- 
tions for the computation of the gravitational potential. 

The initial conditions used in our comparison provide 
a more rigorous test of the performance of our codes than 
those commonly used in the study of cluster formation. 
Typically, initial conditions are set up with roughly equal 
kinetic an d grav itational ener gies, i?kin ~ -E-grav (e.g., 
iBate et al.ir2003HBonnell et al.li2003D . These clouds tend 
to evolve to form several centers of star formation, each 
producing a small group of fragments that feed from the 
gas delivered by the large-scale flows in the cloud. Clouds 
with less kinetic energy, in contrast, have less support 
against global contraction, and so undergo global col- 
lapse. As a result, such clouds tend to produce a sin- 
gle, dense star-forming cen ter that also ha s a higher rate 
of star formation (e.g., Cl ark et al.l 120081 ). As a conse- 
quence, the cloud modeled here enters a more violent 
and chaotic A^-body phase, in which it is more difficult 
to obtain convergence for different numerical methods. 

4.2.1. Resolution Cntenon in FLASH 

In order to satisfy the resolution criterion in simula- 
tions including self-gravity, the gas density on the grid 
must not exceed a critical density Pres i n regions of grav- 
itational collapse ()Truelove et al.|[l997l ). This density is 
related to the smallest resolvable Jeans length Aj on the 
highest level of the AMR hierarchy. The density threshold 
Pics is obtained by solving equation ([2|) for p, and using 
the Jeans length and sound speed of the gas resolved with 
AMR on the highest level of refinement: 



PvGS 



4Gr 



2 

acc 



(32) 



Since the Jeans length should be resolved with at le ast 4 
grid cells in AMR simulations (|Truelove et al . 1 119971 ) . the 
sink particle accretion radius race must not be smaller 
than 2 grid cells. The sink particle accretion radius is 
thus determined by the smallest linear cell size Ax that 
can be resolved with pure AMR. We set race — 2.5Ax 
to satisfy the Truelove criterion on the highest level of 
refinement. The Jeans length Aj = 2racc is thus re- 
solved with 5 grid cells on the top of the AMR hierar- 
chy. However, on all AMR levels smaller than the max- 
imum level of refinement we resolve the Jeans length 
with at least 12 grid cells to follow features on their way 
to runaway collapse more accurately. We additionally 
apply the shock refi nement criterion provided in flash 
(jFrvxell et al.|[2000l), and our additional refinement on 
sink particles (cf. 

For the comparison test presented here, we used a fixed 
base grid with 128^ grid cells plus two levels of AMR 



resulting in an effective resolution of 512^ grid cells. 

4.2.2. Resolution Criterion m SPH 

IBate fc Burkerd (|1997[ ) have shown that at least two 
times the number of SPH neighbors, 2A^noigh particles are 
necessary for resolving the minimum Jeans mass to avoid 
artificial fragmentation in SPH codes. The Jeans mass is 
defined as 



15/2/^ ^bT 



6 \ Gfimu 



3/2 



-1/2 



(33) 

where fee, M ^-nd toh are the Boltzmann constant, the 
mean molecular weight and the mass of a hydrogen atom, 
respectively. Assuming a number of SPH neighbors of 
roughly A'ncigh — 50, we can estimate the mass resolu- 
tion limit for an SPH calculation as Mies — Afj(Pros)/100. 
The corresponding sink particle accretion radius is race — 

Aj(pres)/2. 

Since we resolve the Jeans length with 5 grid cells in 
the FLASH run, we chose to use a similar number of SPH 
smoothing lengths to resolve the Jeans length in the SPH 
run. Therefore, we used a total number of 2 million 
SPH particles for our comparison, which roughly corre- 
sponds to resolving the Jeans length with 5 SPH smooth- 
ing lengths. Furthermore, the initial velocity field was 
constructed on a grid with 128^^ ~ 2 x 10® grid cells. 
Thus, to accurately sample the initial velocity field with 
SPH particles, it was reasonable to use a similar number 
of resolution elements in the SPH code as grid cells for 
the initial velocity field. 

4.3. Results of the Sink Particle Code Comparison 
4.3.1. Column Density Distributions 

FigurefTSlshows the column density distributions of the 
FLASH (left) and the SPH run (right) at t = 0.0 iff (top) 
and t = 0.5 iff (bottom). The initial supersonic turbulent 
velocity field has generated a complex network of shocks 
at t = 0.5 tg. Some of the shock-compressed gas becomes 
gravitationally unstable and goes into freefall collapse. 
This happens at roughly t = 0.8 iff as shown in the top 
panels of Figure [14] (left and middle columns for flash 
and SPH, respectively). Five sink particles have formed 
in both the flash and the SPH run, containing a total 
mass of 0.9 Mq and 0.7 Mq, respectively. They form at 
roughly the same locations in the gas distribution with 
slight differences in the exact positions. These differences 
are mainly due to the slightly different hydrodynamic 
evolution of the gas prior to sink particle formation in 
the FLASH and SPH runs. There is striking agreement 
between the flash and SPH runs at t = 0.9 tg, indicated 
by the number of sink particles (15 versus 17) and their 
accreted mass (4.9 M0). 

At i = 0.9 tff, the gas distribution has already started 
to collapse globally due to the decay of the initial turbu- 
lence. At i = 1.0 iff, 32 sink particles containing a total of 
14.5 Mq have formed in the flash run, while the SPH run 
has produced 48 sink particles with 20.0 A/©. Thus, the 
FLASH and SPH runs start to diverge at t = l.Otg. The 
star cluster produced in the SPH run is slightly denser 
than in the flash run. This trend continues to the last 
snapshot at t = 1.1 tff when the FLASH run has produced 
49 sink particles with 26.3 Mq, while the SPH run has pro- 
duced 61 with 36.7 Mq. The star cluster is clearly denser 
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in the SPH run at this last time. There are still new sink 
particles forming in the outskirts of the cluster. They 
from at very similar locations in both the FLASH and the 
SPH runs. For example, a.t t — l.lts the rightmost sink 
particle has formed at almost exactly the same location 
in both codes with respect to the position of the dense 
core forming there. However, this particular sink particle 
also forms slightly closer to the center of the star cluster, 
reflecting the faster global collapse in the SPH run. This 
is discussed further in the following section. 

4.3.2. Radial Profiles 

To quantify the differences in the flash and SPH runs 
we show radial profiles of the mass, density and radial 
velocity in Figure [H] for times t = 0.0, 0.5, 0.8, 0.9, 1.0 
and 1.1 iff (from top to bottom). The initial conditions 
are compared in the three top panels. The density is 
constant inside the initial gas sphere with a radius of 
i? = 5.0 X 10" cm = 3.3 x 10"* AU. Thus, the integrated 
mass grows with the radius r proportional to inside 
the sphere until it reaches M = 100 Mq at r = i?. In 
the FLASH run, the mass grows slightly further for r > 
R up to about 103 Mq, because of the remaining low- 
density gas outside of the sphere as discussed in § 14.21 
In contrast, the density is zero for r > R in the SPH run. 
Initially, the radial velocity is roughly zero, because the 
turbulent velocity field was constructed to be isotropic 
with zero mean. 

At t = 0.5 iff the radial velocity indicates the onset of 
global infall with slightly supersonic speeds of 0.5kms~^ 
at r ~ 1 - 2 X 10"^ AU. Both the flash and the SPH 
run exhibit similar radial profiles at that time. When 
the first sink particles have formed {t = 0.8 te) the radial 
density profiles roughly follow p{r) cx r^^/"^ (see e.g., 
IT977I) for 2 X 10^ < r/AU < 2 x lO"*. At t = 0.9 iff the 
FLASH and SPH run have produced 15 and 17 sink par- 
ticles, respectively, which is seen as density peaks in the 
radial density profiles. The infall velocities have reached 
Mach 5 at that time, and they are roughly consistent be- 
tween the two codes. However, a.t t = 1.0 iff the SPH run 
shows slightly larger infall speeds than the flash run for 
2 X 10"^ < r/AU < 2 X lO'*. This becomes more promi- 
nent at t — 1.1 iff, reaching down to the cluster center, 
when the infall speeds in the SPH run are almost twice as 
large as in the flash run. The star cluster produced in 
the SPH run is denser, which is shown by the mass and 
density profiles at t — 1.1 iff. This is consistent with the 
visual inspection of column density images discussed in 
the previous section. Thus, the global collapse of the gas 
cloud proceeds slightly faster in the SPH run compared 
to the FLASH run. 

There are three explanations for this behavior, proba- 
bly all contributing to the faster collapse speeds in the 
SPH run. First, the gas mass outside the sphere may con- 
tribute to slow down the collapse slightly in the FLASH 
run. Inflow boundary conditions were used for the hy- 
drodynamics in the flash run. Thus, low-density gas 
is streaming in from the boundaries during the global 
collapse, which further increases the filamentary excess 
mass outside of the cluster from 3 Mq to 5 Mq. Secondly, 
both the FLASH a nd SPH codes have di fferent dissipa- 
tion mechanisms. iKitsionas et al.l ()2009l ) find that SPH 
codes tend to dissipate small-scale, random, turbulent 



motions slightly faster than grid codes. Thus, our SPH 
run may loose turbulent support faster than our flash 
run, which shifts the onset of global collapse to earlier 
times in the SPH run. Finally, one must consider that 
the low-de nsity gas is better resolved in the grid calcula- 
tion (e.g.. lPrice fc Federrath ,2010;) . This becomes even 
stronger at later times when the number of SPH particles 
decreases due to accretion, and the number of grid cells 
in the FLASH run increases due to the Jeans refinement 
and the formation of sink particles (cf. I2.2.2p . There- 
fore, at late times, the numerical resolution of the flash 
run becomes increasingly superior to the SPH run. For 
the star cluster simulation presented here, lower numer- 
ical resolution could lead to somewhat smaller collapse 
timescales, because small-scale turbulent pressure sup- 
port is weaker at lower resolution. We conclude that all 
these three effects may contribute to the slightly faster 
global collapse in the SPH run compared to the flash 
run. 

4.3.3. Mass Distributions of Fragments 

Figure [in] shows a comparison of the mass distributions 
of sink particles obtained in the flash and SPH runs at 
i = 0.8, 0.9, 1.0 and 1.1 iff (from top to bottom). The left 
panels show the histograms of sink particle mass, while 
the right panels show the cumulative mass distributions 
from which the Kolmogorov-Smirnov probability (p-KS) 
was determined with a KS test. The KS test provides an 
estimate of the probability that the mass distributions 
obtained in the flash and in the SPH run were drawn 
from the same basic distribution, and it thus provides an 
estimate of their similarity. 

Figure [12] shows that at i = 0.8 iff the KS probability 
is about 70%, which means that the mass distributions 
obtained in the flash and SPH run are very likely drawn 
from the same distribution. The similarity in the num- 
ber of particles formed and in the mass accreted onto sink 
particles taken together with the similarity in the loca- 
tion of sink particles at that time is encouraging. This is 
an important result, because it shows that two indepen- 
dent implementations of sink particle formation in two 
fundamentally different hydrodynamic codes give almost 
identical results for the number, mass and location of 
formed fragments. 

At later times, the KS probability drops to 55%, 40% 
and 8% at i = 0.9, 1.0 and 1.1 iff, respectively. How- 
ever, the general shape of the mass distributions obtained 
in the flash and SPH runs is similar for all times (al- 
though the relatively small number of sink particles does 
not allow for a detailed comparison of the slope at the 
high- mass end of the distributions). The peaks of the 
mass distributions agree well for i < 1.0 iff. However, 
at i — 1.1 iff the peak of the distribution is shifted to 
slightly higher masses in the SPH run as a result of the 
faster global collapse in the SPH run (cf . § 14.3.21) . The 
high- mass ends of the flash and SPH mas s distri butions 
are roughly consistent with the ISalpeteij (jl955f) power 
law, N cx M~^-^^, as also found in previous SPH studies. 
However, our statistical samples are not large enough to 
draw further quantitative conclusions on the slope of the 
high-mass end. It should also be noted that the low- 
mass end is probably affected by our assumption of an 
isothermal equation of state. Our models did not take 
into account the effects of radiation transfer, which is ex- 
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pected to s uppress fragmentation (jKrumliolz et al.|[2007l : 
lBatdl2009D . However, the fact that our sink particle ra- 
dius is about 500 AU in this simulation does not allow 
any fragmentation below this scale. Thus, we are more 
likely underestimating the amount of fragmentation, be- 
cause we cannot follow possible fragmentation on scales 
smaller than about 500 AU. 

The fact that the KS probability becomes smaller at 
later times is a result of the faster global collapse in 
the SPH run. To quantify this further, we compare the 
mass distributions at the times when the total mass of 
sink particles is roughly the same in the flash and SPH 
runs. Figure [T7] shows the mass distribution of flash 
at t = 1.1 iff' together with the mass distribution of SPH 
a,t t = 1.04 iff (about 5% earlier) when the sink parti- 
cle formation efficiency is about 26% in both runs. The 
mass distributions agree very well for similar sink for- 
mation efficiencies, which is indicated by a KS proba- 
bility of 55%. Moreover, the number of fragments is 
almost identical in the flash and SPH run (49 versus 
50). There is striking agreement between the sink par- 
ticle properties obtained in our flash implementation 
and in the Bate ct al. (1995) implementation for SPH al- 
beit the highly complex and chaotic nature of the prob- 
lem modeled in this code comparison. The fact that two 
fundamentally different numerical schemes with comple- 
mentary strengths and weaknesses give the same overall 
result is very encouraging. It strengthens our confidence 
in simulations of the turbulent collapse of clouds with 
AMR on the one hand and SPH on the other hand. 

4.3.4. Computational Efficiency of FLASH and SPH 

We briefly mention the computational efficiency of 
our FLASH and SPH runs. Both codes were run in 
a mode of parallel computation on the same super- 
computer (HLRB-II: SGI Ahix 4700) at the Leibniz- 
Rechenzentrum Garchin^. The flash run consumed 
about 10,300 CPU hours and was run on 128 CPUs. 
The SPH run consumed roughly 2,400 CPU hours and 
was run on 16 CPUs, which is 4-5 times faster than the 
FLASH run. The SPH code automatically increases reso- 
lution in high-density regions, without the necessity to 
increase the number of resolution elements there, be- 
cause the SPH particles move in the Lagrangian frame. 
It should be noted however that due to accretion of SPH 
particles onto sink particles, the total number of SPH 
particles decreases from the initial 2 million to about 
1.3 million at t = 1.1 iff. In contrast, the Eulerian code 
FLASH refined the grid adaptively in high-density regions. 
The number of resolution elements (grid cells) thus in- 
creases from about 32 million at t = to 56 million at 
t = l.ltff, which is a factor of about 30 more resolution 
elements on average in the flash run compared to the 
SPH run. The computational cost per resolution element 
was thus a factor of 6-8 smaller for flash than for SPH. 
This is roughly consistent with the difference in the aver- 
age computational efficiency per resolution element mea- 
sured for typical grid-based and particle-based hydrody- 
namic codes in t he case of non-selfgravitating, supersonic 
turbu lence fe.g. JKitsionas et "aLll2009HPrice fc FederrathI 
[2OT0h . 

http://www.lrz-muenchen.de 



5. CONCLUSIONS 

We introduced an implementation of accreting sink 
particles for the adaptive mesh refinement code flash. 
Sink particles are used as a subgrid model to cut-off local 
runaway collapse in a controlled way to avoid artificial 
fragmentation and to avoid prohibitive small timesteps 
in numerical simulations involving multiple collapse re- 
gions. Sink particles interact gravitationally with the gas 
and with one another, and can accrete bound gas. The 
gas must pass a series of checks prior to sink particle 
formation as explained in detail in § 12.21 In particular, 
in our grid-based implementation the gas considered for 
sink particle creation is not only required to exceed a 
given density threshold, but this gas must be gravita- 
tionally bound and collapsing at the same ti me, similar 
to the original sink particle implementation of iBate et al.l 
(|1995D for SPH. Figure [m (right panels) shows that spu- 
rious sink particles would be created in the star cluster 
formation run discussed in § |4l if a sole density thresh- 
old was used to determine sink particle creation. The 
root cause of spurious sink particle creation is that the 
gas density can exceed the given density threshold in 
shocks that do not necessarily trigger gravitational col- 
lapse. Figure [T4l fright versus left panels) shows that the 
number of fragments would be overestimated by more 
than one order of magnitude, and that the star forma- 
tion efficiency (total gas mass accreted by sink parti- 
cles) would also be overestimated by 133%, 51%, 19% 
and 14% at t/tg = 0.8, 0.9, 1.0 and 1.1, respectively, if 
the additional checks described in § 12.21 are switched off. 
Some sort of sink particle merging could be used to re- 
duce the overestimated number of sink particles, but the 
mass accreted onto sink particles would still be overesti- 
mated if a sole density threshold was used for sink par- 
ticle creation at densities below the opacity limit (about 
10-i4gcm-3). 

We performed a series of tests of our sink particle im- 
plementation in § [21 including the time integration of 
circular and elliptical orbits, the collapse of a Bonnor- 
Ebert sphere, and a rotating cloud core fragmentation 
test. The sink particle accretion rate in the collapse of 
a singular isothermal sphere showed excellent agreement 
with the theoretical predictions of the Shu (1977) model. 
Those tests showed that the dynamical creation, gravi- 
tational interaction, motion, timestepping and accretion 
of sink particles works in a way that enables us to ana- 
lyze the trajectories, accretion rates and mass distribu- 
tions of collapsing fragments quantitatively and reliably 
in follow-up studies that will also include mhd effects. 

A comparison calcula tion of star clust er formation be- 
tween the SPH code by iBate et al.l (|l995f ) and our grid- 
based flash implementation showed encouraging agree- 
ment of gas and sink particle properties obtained in the 
SPH and flash runs (cf. §[4]). Using column density im- 
ages, we demonstrated that sink particles are created at 
roughly the same locations and times in both codes. Ra- 
dial velocity profiles revealed that the global collapse in 
the SPH run was about 5% faster than in the FLASH run, 
probably due to the slightly higher nu merical d issipation 
in SP H codes (Kitsionas ct al. 2009; IPrice fc Fcdcrrat"hl 
12010'). After correction for this time lag, we showed that 
the number of sink particles and their mass distributions 
in the flash and SPH runs are in very good agreement. 
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The agreement of our flash and SPH runs is an encourag- 
ing result that strengthens our confidence in numerical 
simulations of gravoturbulent cloud fragmentation and 
collapse with AMR and SPH. 
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APPENDIX 
GRAVITATIONAL SPLINE SOFTENING 



Cubic spline softening fsee lMonaghan fc Lattanziolll985l : lPrice fc MonaglL"anll2007| ) is used to soften the gravitational 
accelerations g(r) between sink particles and gas and among sink particles: 
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for i < ^ < 1 

2 — r^alt 



for r > Tsoft , 



(Al) 



where r = |r| and Tgoft is the softening radius. See Figure [T] for a graphical representation of this type of softening 
compared to Plummer softening and compared to Newton's acceleration of a point mass. 

CONSERVATION LAWS DURING ACCRETION OF GAS 

We briefly discuss the conservation laws during accretion of gas onto sink particles. The index i can be used to 
denote both the multiple gas portions to be accreted and the sink particle onto which these gas portions are being 
accreted within a single accretion step. Thus, mass and linear momentum conservation during accretion or merging 
of sink particles are given by the following equations: 



i 



(Bl) 

(B2) 



The last equation determines the center of mass velocity Vcm to conserve linear momentum during an accretion or 
merging step. However, if these two conservation laws apply, then angular momentum conservation is generally broken 
in the process of accretion. To see that, consider the angular momentum before the accretion process 



TO,; r,; X Vi 



The angular momentum after accretion is 



L' = MR X Vp 



(B3) 



(B4) 



To fulfill angular momentum conservation L = L', the following equation needs to be solved for the position R of the 
accreted gas plus sink particle: 

1 



M 



(B5) 
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However, this equality cannot be fulfilled in the general case, because L and Vcm are generally not orthogonal. Thus, 
it is impossible to construct a position vector R, such that angular momentum conservation holds. Furthermore, if R 
is given by the center of mass 

Rem = ]^ X! 

i 

angular momentum conservation is also typically broken, because eq. (|B5|) does not generally hold true. This is the 
case for any implementation of accreting sink particles. The only way to res tore angular momentum conservation is 
to introduce an intr insic angular momentum (spin) for each sink particle (see lBate et al.l [19951 : iKrumholz et al.ll2005 
iJappsen et £01120051 ). This spin compensates the deviation from global angular momentum conservation caused by 
accretion or merging of sink particles. It can also be used to determine the axis along which a bipolar outflow or jet 
would be launched in a subgrid model of stellar feedback. 
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^ 1 



0.0 0.5 




Fig. 1. — Plummer softening versus spline softening; The sink particle implementation presented here uses spline softening (Appendix IXl 
for the gravitational accelerations between sink particles and grid cell centers, and between adjacent sink particles. 




Fig. 2. — Trajectory of two point masses on circular orbits around their common center of mass after 1000 orbits of time integration. 
The leapfrog integrator is symplectic for constant timesteps, and thus the circular orbits remain numerically stable. 
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Fig. 3. — Trajectories of two point masses on elliptical orbits around their common center of mass after 10 orbits. The symplcctic property 
of the leapfrog integrator is broken due to variable timesteps in this type of problem. However, 10 orbits are reproduced quite accurately 
with subcycling turned ON, while for subcycling turned OFF (bottom), large errors occur already for single orbits, and a significaiit amount 
of artificial angular momentum is introduced. 
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Fig. 4. — Gravitational softening and subcycling test: The trajectories of two point masses starting with zero velocity along the j;-axis, 
but with an initial velocity along the a;-axis are shown. The two masses oscillate around their common center of mass, and draw a sine 
function. The boundary conditions are set to outflow for the y-direction and to periodic for the a;-direction. The frequency of the oscillation 
in y and the initial velocity in x were chosen such that the two particles leaving on one side of the domain will exactly connect to their 
previous trajectories. Each sweep through x contains 20 oscillations and a total number of 1000 oscillations is shown (top). When - 
for the same setup - subcycling is turned OFF (bottom) the sink particles get artificial accelerations during their close encounters, and 
energy conservation is broken. Already after roughly 68 oscillations, they have artificially gained enough energy to leave the domain at 
X « -0.1 AU. 
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Fig. 5. — Three sink particles orbiting a static spherical gas density profile with p oc r~'^ . Note that the hydrodynamical evolution was 
switched-ofF deliberately for this test, while the gravity solver and the AMR scheme remained fully operational. The left image shows the 
initial conditions and the right image shows the sink particle positions after the innermost particle has completed roughly one full orbit. 
Adaptive mesh refinement is used to keep the sink particles at the highest level of refinement allowed in this setup. The particle on the 
middle orbit has completed half an orbit, while the outermost particle has completed 1/3 of an orbit at the time shown in the right-hand 
image. Figure [6] shows the full particle trajectories after one, ten and 100 orbits of the innermost sink particle. 
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Fig. 6. — Sink particle trajectories of the setup shown in Figure [5] after one (left), ten (middle) and 100 (right) completed orbits of the 
innermost particle. The time in the left plot is equivalent to the image shown in Figure [5] (right). Ten orbits (middle) are quite accurately 
followed. After 100 orbits of the innermost sink particle (right), the intermediate particle has finished 50 orbits, while the outermost sink 
particle has finished about 33 orbits. The outermost sink particle shows some deviation from a perfect spherical orbit, seen as a slight 
broadening of its circular trajectory, which is due to the accumulation of small time-integration errors, and due to deviations from spherical 
symmetry toward the edges of the Cartesian domain. 
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Fig. 7. — Time evolution of the mass accretion rate, dm/dr, and of the accreted mass, m, of the central sink particle in the case of a 
collapsing Bonnor-Ebert sphere. The quantities are given in dimensionless units (see text). The mass accretion rate increases rapidly until 
about 12% of the gas mass is accreted onto the sink particle (the total mass of the Bonnor-Ebert sphere is m = 17.3). Here, the density 
threshold for sink particle accretion was set to 10^ po, which corresponds to an effective accretion radius of J = 0.1. The time Tsink = 
corresponds to the time when the sink particle starts to accrete. 




0.25 



Fig. 



Same as Fig. [7] but in the case of prcs = 10* po, corresponding to an effective accretion radius of J = 10 ^, ten times smaller 



than in the previous test. The peak accretion rate is thus reached ten times faster than in Fig. [7] 
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time [yr] 

Fig. 9. — Mass accretion of the central sink particle during the collapse of a singular isothermal gas sphere. The mass accretion rate onto 
the sink particle is roughly constant in time with a value of M 1.5 X 10~'^Mq yr~^, until the sphere runs out of gas at t 2 X 10'* yr. 
The collapse is highly supersonic due to the large instability parameter A = iirG p{R) /c^ = 29.3. Therefore, the mass accretion rate 
reaches values much larger than the expansion-wave co llapse solu tion bv lShul 1119771) that would give moc^/G = 1.06 X 10~^ Mp , yr~* for 
the asymptotic solution: A = 2 and mo = 0.975 (see ISliuiri977l . Tab. 1). For our case, A = 29.3 and mo 133, IShul II1977I ) predicts 
an accretion rate of M fa 1.45 X 10~* Mq yr~*, which is in excellent agreement with our numerical estimate of the central sink particle 
accretion rate. 
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Fig. 10. — Time evolution of the radial density profile (top panel), infall velocity (middle panel) and mass flux through radial shells 
(bottom panel) for the collapse of a singular isothermal sphere. The vertical dotted line marks the accretion radius of the sink particle. 
The density profll e evolves from the initial profile to r~^'^ when the mass of the sink particle starts to dominate the gravitational 
potential (see e.g. IShu l 119771 : fOgino et al.lll999h . As expected for a highly unstable sphere the collapse proceeds with highly supersonic 
speeds (note that Cs = 0.166kms~^). The mass fiux of gas through spherical shells of radius r is constant in time. It is exactly the mass 
accretion rate of the central sink particle, M ^ 1.5 X lO~^M0yr~^ (cf. Fig.|9]l. 
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Fig. 11. — Coefficient mo for the accretion rate in equation II29I I as a function of th e insta bility parameter A = A-kG p(R) ^f ^/Cg for the 
collapse of a singular isothermal sphere. The diamonds show tabulated values bv lShul l)1977l . Tab. 1). For A = 29.3 as in the collapse test 
shown in Figures l9l and [TOl the accretion rate is more than two orders of magnitude higher than cf/G, in excellent agreement with the 
predicted accretion rate bv lShul 1)19771 ). 
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Fig. 12. — Column density time sequence of the inner 400 AU of a rotating cloud core (see initial conditions in ^ 13.51 1. the so called 
'BB test' ijBoss fc Bodenheimeii 0-9791 : [Bate &: Burkertl 119971) . Two sink particles form at the location of the two cloud fragments at 
t = 4.31 X 10* yr. At t = 4.39 x 10"* yr a bar forms that connects the two main fragments, in which a third particle is dynamically created 
at t = 4.41 X 10** yr right in the center. The two main fragments move toward the central object, and at t = 4.52 X 10'* yr two of the three 
sink particles merged to a single particle close to the center of the rotating cloud core, because the optional sink particle merging was used 
in this test. The freefall time of the initial gas distribution is = 3.41 X lO'* yr. Sink particles are drawn as black circles with a radius of 
^succ = 39 AU in a true-to-scale representation. 
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Fig. 13. — Comparison of the column density distributions obtained with the AMR code FLASH (left) and with the SPH code bv lBate et al.l 
(|1995l) (right) at t = 0.0 iff (top) and t = 0.5 tg (bottom) prior to the formation of a stellar cluster. The initial supersonic, turbulent 
velocity field creates a complex network of filaments. Some of these shocked regions become gravitationally unst able and collapse, i.e., sink 
particles are created, while other dense regions are transient features, not leading to local collapse (see Fig. I14|l . 
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Fig. 14. — Figure [131 continued. White filled circles show the projected position of sink particles obtained with FLASH (left panels) and 
SPH (middle panels) as a function of the global freefall time {t/tg = 0.8, 0.9, 1.0, 1.1, from top to bottom). Sink particles form at roughly 
the same positions in the FLASH and SPH runs. The number and mass of sink particles is in excellent agreement for t < Itg. At later times 
(t > Itff), the stellar cluster produced in the SPH run collapses on slightly smaller timescales than the FLASH cluster, and is more centrally 
condensed at t = 1.1 ig. The right panels show the effect of only using a density threshold for sink particle creation with FLASH, compared 
to the default FLASH run (left panels), for which all sink particle creation checks (see Section 12.211 were used. Without additional checks 
for converging flows, gravitational potential minimum, bound state, and Jeans instability, a signiflcant number of spurious sink particles 
are created, which do not represent collapsing objects. Most of the sink particles in the right hand panels were formed in transient density 
peaks caused by shocks that did not accumulate enough mass to create gravitationally bound and collapsing objects. The much larger 
number of particles created in that case may be reduced by some sort of merging algorithm, but the accreted mass in sink particles would 
still be overestimated, in particular at early times (by 133%, 51%, 19% and 14% compared to the default FLASH run at t/tff = 0.8, 0.9, 1.0 
and 1.1, respectively), if a sole density threshold was used for sink particle creation at densities below the opacity limit. In contrast, the 
sink particles in the left hand plots (default FLASH runs) represent gravitationally bound and collapsing objects, as in the SPH run (middle 
panels). 
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Fig. 15. — Comparison of radial profiles (cumulative mass, density, and radial velocity from left to right), obtained in the FLASH and SPH 
runs as a function of global freefall time (t/tfi = 0.0, 0.5, 0.8, 0.9, 1.0, 1.1, from top to bottom). 
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Fig. 16. — Comparison of the sink particle mass distributions (left) and cumulative mass distributions (right), obtained with FLASH and 
SPH as a function of global freefall time (t/tg = 0.8, 0.9, 1.0, 1.1, from top to bottom). 
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Fig. 17. — Same as the bottom panels of Figurc FlGl but the mass function of FLASH at t = l.lOfff is compared to the mass function of SPH 
at t = 1.04 iff, when the accreted gas mass is roughly the same for both runs (M* ~ 26 Mq). The number of fragments (sink particles) is 
in excellent agreement between the two codes. The mass distributions are also very similar as shown by a Kolmogorov-Smirnov probability 
of 55%, which indicates that the two samples are likely drawn from the same fundamental distribution. 



